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be achieved with different parameters. The obtained results may raise the possibility of relative 
experiments and potential applications. 
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I. INTRODUCTION 

Solitons describe a class of fascinating nonlinear wave 
propagation phenomena appearing as a result of balance 
between nonlinearity and dispersion or diffraction prop- 
erties of the medium under nonlinear excitations, which 
leads to undistorted propagation over extended distance 
[l[ . One of the most important physically relevant real- 
izations of solitons is provided by the matter- wave soli- 
tons in Bose-Einstein condensed atomic gas Q. Based 
on the successful experimental realization and theoret- 
ical analysis of Bose-Einstein condensations (BECs) in 
weakly interacting atomic gases 0], matter- wave dark 
solitons Q, vortices bright solitons @, gap soli- 
tons 0], and soliton chains Q have been observed and 
studied. These studies have stimulated a large amount of 
research activities, which enable the extension of linear 
atom optics to nonlinear atom optics 

The realization of higher-dimensional matter-wave 
solitons in BECs is still a challengeable topic be- 
cause those solutions are usually unstable for (2+l)-D 
and (3+l)-D constant-coefficient nonlinear Schrddingcr 
(NLS) equation due to the weak and strong collapse |9j. 
However, different situations are observed in BECs with 
temporally or spatially modulated parameters. Alter- 
ation of atomic scattering length achieved by Feshbach 
resonance [10( has been used to dynamically stabilize 
higher-dimensional bright solitons 11| while periodic ex- 
ternal potentials achieved by optical lattice has been used 
to g enerate and control higher-dimensional gap solitons 
ID periodic wave solutions are also predicted in 
BECs with time-space varying parameters 13J. More- 
over, the bright solitons [14| and periodic wave solu- 
tions }15| were obtained in spinor BECs governed by a 
system of three coupled mean-field equations. 

In this work, we present a detailed study on dynamics 
of analytical 3D bright matter-wave single solitons and 



soliton-pairs in BECs with time-space modulation. We 
note that 3D periodic wave solutions have been studied 
in the generalized NLS equation very recently 3, 17 1. 
However, the authors did not study the soliton pair solu- 
tions and their interaction properties. By using the sim- 
ilarity [l3, 17, IH and bilinear transformations we 
can achieve different shapes of bright solitons and fasci- 
nating interactions between two solitons. In addition, the 
experimental possibilities for observability are discussed 
and the stability of solitons is illustrated numerically. 

The paper is organized as follows. In the next section, 
the model under study is introduced. In Sec. Ill, the 
methods for solving the model equation are introduced. 
A relationship between the model and a practical system 
is established. In Sec. IV, we give the expressions of the 
bright solitons and the soliton-pairs. The interactions 
between two solitons are also investigated. In the last 
section, the case of the dark solitons is discussed and the 
outcomes are summarized. 



II. THE GP MODEL 

The dynamics of a weakly interacting Bose gas at zero 
temperature is well described by the (3+l)-D GP model 
with time-space modulation Q 



lh -m = 



-—V 2 +V cxt (t,r) + G(t)\y\ 
2m 



tf + ir(t)¥, (i) 



where V = (d x ,d y ,d z ), r = (x,y,z), * = *(f, r) de- 
notes the order parameter with N = J \ ^\ 2 dr being the 
number of atoms in the condensate, G(t) = ATrh 2 a s (t)/m 
is the interaction function with a s (t) being the s-wave 
scattering length modulated by a Feshbach resonance, 
and T(t) is the gain/loss term, which is phenomeno- 
logically incorporated to account for the interaction of 
atomic or thermal clouds. We note that the dissipative 



2 



dynamics originating from the interaction between the 
radial and axial degrees of freedom has also been stud- 
ied recently [2(|. Here the potential is chosen as a har- 
monic trap V cxt (t,r) = ( m /2)(r — e(t))uj 2 (t)(r — e(t)) 
with uj(t) — dia,g(oj x (t),Wy(t),oj z (t)) being a diagonal 
matrix of the trap frequencies in three directions and 
e(t) =(ei(t), e2(i), e 3 (t)) corresponding to its center. 

Using the suitably scales and variables: r = a z r' , t = 
Q-H', e(t) = a z e'(t), = y/N/afy, a z = [fi/K)] 1 ' 2 , 
and oj z — f u) z (t)dt, we arrive at the dimensionless GP 
equation in the (3+l)-D space after dropping the primes 



i. 



v(t,r) + fl (t)H s 



^ + i 7 (t)V, (2) 



where g(t) = 4wN 2 a s (t)/aj, *y(t) = T(t)/(huj z ), and 



v(t,r) = -(r-e(t))a 2 (t)(r-e(t)) 



(3) 



with a(t) = diag(a 1 {t),a 2 (t),a 3 (t)) = wj 1 ^). Eq. © 
is associated with 5£/8ip* = in which the Lagrangian 
density can be written as 



-2[«(t,r)+fy(t)]|Vf. 



(4) 



Here for the harmonic trapping potential v(t,r) given 
by Eq. ([3]) , after some algebra it follows from system ([7]) 
that the similarity variables can be expressed as 



£(t,r)=/?(t)-r+ / P(s) ■ a(s)ds, 



(8) 



r(i)= / \mrda, 
Jo 



where (3(t) — (fii(t), (3 2 (t), (3 3 (t)) denotes the vector 
of the inverse spatial widths of the localized solu- 
tions along x, y, z directions, respectively, and cr(t) — 
(ai(t),a 2 (t),a 3 (t)) with aj = (3j J ejO^fiJ 1 dt relating 
to the velocity of the solitons. Moreover the nontrivial 
phase has the quadratic form 



<p(t,T) = -±TA(t)T + v-T-^J ( 



ea 2 e 



W\ 2 )dt, (9) 



where A(t) = diag(/3i//3i, fc/fa, k/fa). The additional 
relations between ctj (t) and (3j (t) = 1/ Vj (t) result in the 
Mathieu equations 



i> j (t)+a 2 Jt)v j (t) = 0, (j = 1,2,3). 



(10) 



Finally, the function p(t) modulating the amplitude of 
solution ip and nonlinearity g(t) can be also found by 



III. SIMILARITY SOLUTIONS 

Here we focus on the spatially localized bright solitons 
and soliton-paris for which lini| r |_ ! . 00 ?p(t, r) = 0. Our 
first objective is to reduce Eq. ([2]) to the tractable NLS 
equation 



.9$(r,0 ld 2 $(T,£) 



+ Gmr,0\ 2 Hr,0 (5) 



using a proper similarity transformation, where r = r(i) 
and £ = £(i, r) are both the unknown variables, and Q is a 
constant. We explore the attractive nonlinearity, i.e. Q = 
— 1, resulting in the bright multi-soliton solutions. The 
case Q = 1 resulting in the dark multi-soliton solutions 
does not pose new challenges and will be discussed in the 
last section. Using the similarity transformation (l3. 17. 
M 



^(t ) r)= /0 (t)e^ r )$(r(t),£(t ) r)), 



(6) 



and requiring $(r(i), r)) to satisfy Eq. and i[>(t, r) 
to be the solution of Eq. ([2]), we find a set of equations 



V 2 £ = 0, & + V£-V^ = 0, r^|VC| 2 
1 



(7a) 



v(t,r) = ~\Vv\ 2 - Vu g(t) = Gp- 2 \Ve, (7b) 



7W = ^f + Pt/p- 



(7c) 



p(t) = Po [f3 1 (t)f3 2 (t)l3 3 (t)}^ 2 exp 
g(t) = -p- 2 {t)\p{t)\\ 



r y(s)ds 



(11) 



which depend on both /3j(t) and gain/loss coefficient -f(t) 
with po being a non-zero parameter. Note that for the 
given /3j(t), the nonlinearity g(t) must attenuate (grow) 
exponentially in the gain (loss) medium ^(t) > (< 0). 

For the given ttj(t), one can, in principle, obtain cor- 
responding f3j(t) (or cquivalcntly for the given Pj(t) 
one can obtain ctj(t)) based on Eq . (p~0|) . Furthermore, 
the bright 7V-soliton solutions of Eq.© can be ob- 
tained using the bilinear transformation [l9| : $at = 
P {N) (t^)/Q {N) (t,0- Here, pW and satisfy 
(iD T + 1/2£ 2 )PW • QW = and D 2 Q^ ■ QW = 
2|p(.iV)|2 ^^jj ]j t an j being the bilinear operators 
and PW = E^ie 2 ^ 1 ^!^^) and QW = 1 + 

SjLi e23 Q2j(r,0- Thus, by choosing /3j(t) and j(t), we 
can generate r) and g(t) for which the generic bright 
./V-soliton solutions of Eq. ([2]) can be found from Eq. 
([5|) on the basis of Eq. ([6|). We will use this analytical 
result to construct the exact bright iV-soliton solutions 
with many interesting nontrivial features. 

For the convenience of analyzing different dynamical 
regimes described by the given model, we specify the 
magnitude of main physical parameters, which are fea- 
sible in experiments. We consider a condensed sodium 
sample trapped in the state \ZSi/2,F — 1, mjr = — 1), 
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FIG. 1: (color online). Curves of ay(t), g(t), and 7(t) given 
by Eqs. g2j) and (JT3J) vs t for p = 6.0, 70 = m = 0.1, and 
6 3 =n = 1.0 (7 = 1, 2, 3). 



which has the scattering length a s = 2.75 nm 21]. The 
other parameters can be taken as N = 1.2 X 10 6 and 
Cj z = (2ir) x 21 Hz, which leads to a z — 4.55 and 
yjN/al = 1.13 x 10 2 ^m~ 3 / 2 . To make sure the frequen- 
cies Oij{t) and nonlinearity g(t) are bounded for realistic 
cases, we choose (3j(t) and the gain/loss coefficient j(t) 
as the periodic functions 

P(t) = dn(t, m)b, 7 (t) = 7o ca(t, n), 70 G R (12) 

where b = (61, 62,63) is a real constant vector describ- 
ing the inverse of the width of the potential and the 
frequency, to 6 [0, 1) and n € [0, 1] are the mod- 
ules of Jacobi elliptic functions. It is easy to see that 
j(t) — 7o sech(t) > corresponds to the dissipative case 
when n = 1 and 70 > (we will focus on this condition 
next). In practical systems, the modulations of (3j(t), 
g(t) and j(t) depend on the use of the optical lattice and 
Feshbach-resonance techniques, i.e. we can achieve T(t) 
and a s (t) by exerting particular time-dependent optical 
field and magnetic field. 

It follows from Eqs. (|10p and (|12l) that ctj(t) is given 
by 

a 2 (t) = to 2 [(2 - m 2 ) sd 2 (t, to) - nd 2 (t, m)] . (13) 

Figure 1 shows the curves of ctj(t), g(i), and 7 (i) vs t. For 
simplicity, we take bj = 1, i.e. ot\(t) = 02 (t) = Q^it) cor- 
responding to the isotropic potential, and consider that 
the center of the potential locates at the origin (e 3 = 0). 
We can also change bj to get an anisotropic potential and 
use nonzero ej to obtain moving bright solitons as those 
discussed in ID case [13j . 



IV. BRIGHT SOLITONS AND SOLITON PAIRS 

Based on the discussions in the previous section, we 
arrive at the fundamental 3D time- varying bright solitons 




FIG. 2: (color online). Propagations (left) and contour plots 
(right) of density for the bright soliton (|14[) in (t, £ = b ■ r)- 
space for ri = 0.5, si = po = 2.0, and 70 = 0.01. The 
others are the same as Fig.l. (a) The breathing bright soliton 
propagating in a zigzag trace for m = 0.1. In the given sodium 
sample, the maximum density and width of the soliton are 
about 1.28 jitm and 12.13 fim. The period is about 53.1 
ms. (b) The bright soliton chain for m = 0.9. The maximum 
density and width of the soliton are about 1.28 /1m" 3 and 9.1 
fim. The period is about 91.0 ms. 



ipx (t, r) = np(t) sech [n (£(t, r) - sir(i)) -In |2n|] e lf> , (14) 

where 9 = Si£(t, r) + (rf — s\) /2T(t) + tp(t, r) with n, Sj G 
E, and p(t), ip(t, r), ^(t, r), and r(t) are given by Eqs. @, 
© and (n). 

Figure 2 exhibits the dynamics of the time-varying 
bright soliton (IT4l . A breathing behavior is also evi- 
dent, which can be managed by /3j(t), Jj(t), and pit). 
For the case to = 0, we have fij — bj and ev., = 0, in 
which the travelling-wave bright soliton is obtained. In 
experiments, it can be simply realized for the zero lin- 
ear potential. The bright soliton propagates in a zigzag 
trace for to = 0.1 [see Fig. 2(a)]. An important fea- 
ture is that while to —¥ 1 (7^ 1) resulting in the larger 
period of /3j(t) given by Eq . (li"2"l) . the amplitude of the 
soliton close to the corners attenuates rapidly so that 
a soliton chain is generated [see Fig. 2(b)]. In experi- 
ments, it can be realized by taking 

(jJ x — (jjy — (jJ z — 

to[(2 - TO 2 )sd 2 (i, to) - nd 2 (£, to)] 1 / 2 ^. 

The interaction of the bright solitons plays an impor- 
tant role in the study of BECs. Here, we will also study 
the interaction between two bright solitons. The analyt- 
ical 3D time- varying bright soliton pairs read 



^ 2 (t,r)=p(i) e ^ t ' r )p(<,r)/g(t,r), 



(15) 



where P{t,r) and Q(t,r) can be expressed as the series 
of exponential functions of (t, r) 

2 2 

3=1 j,k=l,j^k 

(16) 

Q(t,r) = 1+ A 3ke rh+vl +ne r > 1+r > 2+r >' +r >\ 

3,k=l 

with r)j = fijtit, r) + %p 2 T(t), fij = rj + is, {r^s^Sj G 
R), A ifc = djS^ifij + n* k )~ 2 , \jk = (pk - fij)[Sk^jk(pk + 
Hi)' 1 ~ 5 3 K kk {^ + Mfc)' 1 ] 0' ^ an d O = 

(|5i^2|) _1 |/ii-Ai2| 2 (AiiA 2 2VAi2A2i-A 12 A 2 iVAiiA22). 

The dynamics of the 3D time-varying bright two- 
soliton solutions (|15|) is exhibited in Figure 3. Under 
the different parameters, we exhibit three cases for two 
weak zigzag solitons without interaction [see Fig. 3(a)], 
two strong zigzag solitons with interaction [see Fig. 3(b)], 
and strong-weak zigzag solitons with interaction [see 
Fig. 3(c)]. Notice that similar with the bright solitons 
shown in Fig. 2(b), for the case m 1 (m ^ 1), the 
amplitudes of the soliton-pairs close to the corners will 
almost decrease to zero so that panel (a) will degenerate 
to two parallel soliton chains while panels (b) and (c) 
will degenerate to the ><-shaped soliton chains. The 
experimental realization of the dynamics regimes for the 
two-soliton solutions is similar with that for the bright 
one-soliton solutions. 

We stress that the important feature that distinguishes 
our solutions from the reported in the literature [3rl6j is 
the appearance of the time- and space-dependent func- 
tions in both the phase and the amplitude and which 
strongly affect the form and the behavior of bright soli- 
tons and their interactions. 

In order to check the stability of the time-varying 
bright soliton (|14|) . we make numerical simulations of 
Eq. ([2|) with the initial conditions given by Eq. (fT4|) 
and different values of m. We find that the bright soli- 
tons are very stable for m being small (e.g. m — 0.1) [see 
Fig. 4(a)]. With the increase of m, the bright solitons 
become unstable [see Fig. 4(b)]. This is because large m 
results in stronger oscillations of f3j(t), which affect the 
coefficients of Eq. © and the behavior of the solutions. 

V. DISCUSSIONS AND CONCLUSIONS 
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FIG. 3: (color online). Propagations (left ) and contour plots 
(right) for collisions between bright two solitons (|15[) in (t, £ = 
b ■ r)-space for 70 = 0.01, <5i,2 = po — 1.0, and m — 0.6. The 
others are the same as Fig.l. (a) Two zigzag solitons without 
interaction for n = 1.0, T2 = 1.2 and Si,2 = 0. In the given 
sodium sample, the maximum density and width of the left 
(right) soliton are about 1.8 /im~ 3 and 7.6 /im (1.5 /im~ 3 and 
9.1 jum). (b) Two strong zigzag solitons with interaction for 
ri,2 — si — — S2 = 1.2. The maximum density and width of 
the left (right) soliton are about 6.2 /im -3 and 3.0 /im (17.5 
/im -3 and 2.4 (im) (c) Two strong-weak zigzag solitons with 
interaction for n = 1.2, T2 = 1.5, si = 1.05 and S2 = 0. 
The maximum density and width of the left (right) soliton 
are about 6.2 /zm -3 and 4.2 ^m (13.1 /im -3 and 3.8 ^m). 
The period of the zigzag oscillation is about 26.5 ms in all 
panels. 



For completeness, we consider the repulsive nonlinear- 
ity in Eq. ([5]), i.e. Q = — 1. In this case, the equation 
admits 3D dark soliton solutions in the form with a non- 
trial phase 

ip(t,r) = {iv + kt'Anh[k{£(t,r) - vT(t))}}e- lt * T(t) , (17) 

where is the chemical potential, v = y/fi — k 2 , and k is 
a free parameter satisfying k 2 < p:. 



In summary, we have analytically constructed the 
novel 3D time-varying bright multi-soliton solutions for 
the (3+l)-D GP equation with time-space modulation. 
We focus on the bounded potential, nonlincarity, and 
gain/loss case to analyze the dynamics of the breathing 
and the zigzag propagation trace of the obtained solitons. 
Different shapes of the one-soliton solutions and the fasci- 
nating interactions between soliton-pairs were achieved. 
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FIG. 4: (color online). Numerical simulations of bright soliton 
from Fig. 2 vs £ at t = 0, 10, and 20. The initial conditions 
are given by Eq. (I14p with si = 0.1, m = 0.1 in panel (a) and 
m — 0.9 in panel (b). The other parameters are the same as 
those used in Fig. 2. 



The stability of bright solitons have been checked nu- 
merically. The method we present here can be extended 



to study the higher-dimensional bright soliton solutions 
of other nonlinear systems and their various interaction 
properties. The model ((2]) can also be extended to de- 
scribe 3D nonlinear optical media with varying coeffi- 
cients [l6| after the transformation z t. The results 
we obtained may raise the possibility of relative experi- 
ments and potential applications. 
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